# SETUP 

## Packages 
library(stargazer)
library(tidyverse)
library(conflicted)
library(ggsci)
library(ggpubr)
library(broom.mixed)
library(here)
library(dotwhisker)
library(tictoc)
library(glmmTMB)
library(haven)
library(parallel)


conflict_prefer("select", "dplyr")
conflict_prefer("filter", "dplyr")

# Cores for parallel processing 
cores <- detectCores() - 2

# Load and Clean data -------------------------

## Load Data 

load(here("data/lapop_10_18.Rdata")) # slim data -- recodes only

# DATA PREP 

# Select Variables
vars_close_cong <- c(
  "close.cong", "dem.best", "sat_dem", "winner", "non_voter", "pol_interest",
  "edu", "wealth", "rural", "wealth", "female", "trust.people", "neigh.safe",
  "econ.eval", "dem.best.tj", "pop.speech", "maj", "v2x_libdem", "vdem.polar",
  "v2x_corr", "mass.polr", "country", "year", "year.fac"
)


vars_close_court <- c(
  "close.court", "dem.best", "sat_dem", "winner", "non_voter", "pol_interest",
  "edu", "wealth", "rural", "wealth", "female", "trust.people", "neigh.safe",
  "econ.eval", "dem.best.tj", "pop.speech", "maj", "v2x_libdem", "vdem.polar",
  "v2x_corr", "mass.polr", "country", "year", "year.fac"
)


dichot.vars.cong <- c("close.cong", "rural", "female", "non_voter", "winner", "year")
dichot.vars.court <- c("close.court", "rural", "female", "non_voter", "winner", "year")


# Select vars and make factors
cong_dat <- lapop.cleaned %>%
  mutate(across(all_of(dichot.vars.cong), factor)) %>% 
  zap_formats() %>% 
  zap_labels()


court_dat <- lapop.cleaned %>%
  mutate(across(all_of(dichot.vars.court), factor)) %>% 
  zap_formats() %>% 
  zap_labels()


# Models to estimate ----------------------------


# Vector of countries (iso3c)
countries <- unique(select(cong_dat, "iso3c"))
countries <- as.vector(countries[,1])
cong_dat$iso3c <- as.character(cong_dat$iso3c)


### Models to estimate ###

# Congress 
cong_mod1 <- as.formula("close.cong ~ dem.best + sat_dem + edu + wealth + rural + wealth + female + 
                      trust.people + neigh.safe + econ.eval + winner + non_voter + pop.speech +
                      maj + v2x_libdem + mass.polr + dem.best.tj + (1 | country:year) + 
                      (winner + non_voter | country) + year.fac")

cong_mod2 <- as.formula("close.cong ~ dem.best + sat_dem + edu + wealth + rural + wealth + female + 
                      trust.people + neigh.safe + econ.eval + non_voter + pop.speech * winner +
                      maj + v2x_libdem +  mass.polr + dem.best.tj + (1 | country:year) +
                      (winner + non_voter | country) + year.fac")

cong_mod3 <- as.formula("close.cong ~ dem.best + sat_dem + edu + wealth + rural + wealth + female + 
                      trust.people + neigh.safe + econ.eval +  pop.speech * winner + pop.speech * non_voter + 
                      maj + v2x_libdem +  mass.polr + dem.best.tj + (1 | country:year) + 
                      (winner + non_voter | country) + year.fac")

# Court
court_mod1 <- as.formula("close.court ~ dem.best + sat_dem + edu + wealth + rural + wealth + female +
                        trust.people + neigh.safe + econ.eval + pop.speech + winner + non_voter +
                        maj + v2x_libdem +  mass.polr + dem.best.tj + (1 | country:year) +  
                        (winner + non_voter | country) + year.fac")

court_mod2 <- as.formula("close.court ~ dem.best + sat_dem + edu + wealth + rural + wealth + female +
                        trust.people + neigh.safe + econ.eval + non_voter + pop.speech * winner +
                         maj + v2x_libdem +  mass.polr + dem.best.tj + (1 | country:year) + 
                        (winner + non_voter | country) + year.fac")

court_mod3 <- as.formula("close.court ~ dem.best + sat_dem + edu + wealth + rural + wealth + female +
                        trust.people + neigh.safe + econ.eval + pop.speech * winner + pop.speech * non_voter +
                         maj + v2x_libdem + mass.polr + dem.best.tj + (1 | country:year) +
                        (winner + non_voter | country) + year.fac")


# Loop: Drop country, estimate model, and save results -------------------------

# --------------------------------------------------------------------#
# For each country: 
# 1) Drop country from the data 
# 2) Estimate close cong models (1-3) and close court models (1-3)
# 3) Save to lists created below
# --------------------------------------------------------------------#

cong_models1_it <- list()
cong_models2_it <- list()
cong_models3_it <- list()
court_models1_it <- list()
court_models2_it <- list()
court_models3_it <- list()

tic()

for(i in 1:length(countries)) {
  
  # Drop country i 
  new_cong <- subset(cong_dat, iso3c != countries[i]) 
  new_court <- subset(court_dat, iso3c != countries[i]) 
  
  # Tracker 
  print(paste("dropping country:", countries[i]))
  
  # Close Congress Models
  new_cong_fit1 <- glmmTMB(formula = cong_mod1, data = new_cong,
                         family = binomial(link = "logit"), 
                         control=glmmTMBControl(parallel = cores))
  
  new_cong_fit2 <- glmmTMB(formula = cong_mod2, data = new_cong,
                         family = binomial(link = "logit"),
                         control=glmmTMBControl(parallel = cores)) 
  
  new_cong_fit3 <- glmmTMB(formula = cong_mod3, data = new_cong,
                         family = binomial(link = "logit"), 
                         control=glmmTMBControl(parallel = cores)) 
  
  # Close Court Models
  new_court_fit1 <- glmmTMB(formula = court_mod1, data = new_court,
                          family = binomial(link = "logit"),
                          control=glmmTMBControl(parallel = cores))  
  
  new_court_fit2 <- glmmTMB(formula = court_mod2, data = new_court,
                          family = binomial(link = "logit"), 
                          control=glmmTMBControl(optimizer=optim,
                                                 optArgs=list(method="BFGS"), 
                                                 parallel = cores))
  
  new_court_fit3 <- glmmTMB(formula = court_mod3, data = new_court,
                          family = binomial(link = "logit"),
                          control=glmmTMBControl(optimizer=optim,
                                                 optArgs=list(method="BFGS"), 
                                                 parallel = cores))
  
  
  # Save cong models
  cong_models1_it[[i]] <- new_cong_fit1
  cong_models2_it[[i]] <- new_cong_fit2
  cong_models3_it[[i]] <- new_cong_fit3
  
  # Save court models
  court_models1_it[[i]] <- new_court_fit1
  court_models2_it[[i]] <- new_court_fit2
  court_models3_it[[i]] <- new_court_fit3
  
}

toc()

# Save models for later
# Each of these objects is a list of 18 models
save(cong_models1_it, cong_models2_it, cong_models3_it, file = here("data/cong_models_iterative.Rdata"))
save(court_models1_it, court_models2_it, court_models3_it, file = here("data/court_models_iterative.Rdata"))

# Load saved objects 
load(here("data/cong_models_iterative.Rdata"))
load(here("data/court_models_iterative.Rdata"))

# 3 non-positive definite matrices  

## Check the model convergence ----------------------------

# Loops over models in lists 
for (i in 1:18) { 
  
  # Congress models 
  if(cong_models1_it[[i]][["fit"]][["convergence"]] == 1) {
  print(paste("cong 1:", countries[i]))  }
  
  if(cong_models1_it[[i]][["fit"]][["convergence"]] == 1) {
    print(paste("cong 2:", countries[i]))  }
  
  if(cong_models3_it[[i]][["fit"]][["convergence"]] == 1) {
    print(paste("cong 3:", countries[i]))  }
  
  # Court models
  if(court_models1_it[[i]][["fit"]][["convergence"]] == 1) {
    print(paste("court 1:", countries[i]))  }
  
  if(court_models2_it[[i]][["fit"]][["convergence"]] == 1) {
    print(paste("court 2:", countries[i])) }

  if(court_models3_it[[i]][["fit"]][["convergence"]] == 1) {
    print(paste("court 3", countries[i]))  }
  
}


# Tidy Model Results  --------------------------

cong_dfs <- list()
court_dfs <- list()
vars <- c("pop.speech:winner1", "pop.speech:non_voter1", "pop.speech")

# -----------------------------------------------------------------#
# For each specification:
# - Collapse the list of 18 mods into a single dataframe
# - Label the specification 
# - Keep only the main vars 
#------------------------------------------------------------------#

for (i in 1:3) {
  
  #Close congress 
  cong_dfs[[i]] <- map_dfr(get(paste0("cong_models", i, "_it")), tidy) 
  
  cong_dfs[[i]]$submodel <- i
  
  cong_dfs[[i]] <- filter(cong_dfs[[i]],term %in% vars)

  #Close court   
  court_dfs[[i]] <- map_dfr(get(paste0("court_models", i, "_it")), tidy)
  
  court_dfs[[i]]$submodel <- i + 3
  
  court_dfs[[i]] <- filter(court_dfs[[i]], term %in% vars)
  
}

## Label the models with the country dropped ------------------------

country_dropped_m2 <- c()
country_dropped_m3 <- c()

for (i in 1:18) {
  
  # 2 coefficients in model 2 
  new_value_m2 <- rep(countries[i], 2)
  country_dropped_m2 <- c(country_dropped_m2, new_value_m2)
  
  # 3 coefficients in model 3
  new_value_m3 <- rep(countries[i],3)
  country_dropped_m3 <- c(country_dropped_m3, new_value_m3)
}

# Label congress models 
cong_dfs[[1]]$model <- countries
cong_dfs[[2]]$model <- country_dropped_m2
cong_dfs[[3]]$model <- country_dropped_m3

# Label court models
court_dfs[[1]]$model <- countries
court_dfs[[2]]$model <- country_dropped_m2
court_dfs[[3]]$model <- country_dropped_m3


## Make into dataframes and relabel the coeffs -----------------------------

# Dataframes 
cong_models <- rbind(cong_dfs[[1]], cong_dfs[[2]], cong_dfs[[3]])
court_models <- rbind(court_dfs[[1]], court_dfs[[2]], court_dfs[[3]])

# Label the coefficients 

cong_models <- cong_models %>% 
  mutate(term = case_match(term, 
                          "pop.speech" ~ "Populist Discourse", 
                          "pop.speech:winner1" ~ "Pop. Discourse \nx Pres. Supporter",
                          "pop.speech:non_voter1" ~ "Pop. Discourse \nx Non-Voter"))


court_models <- court_models %>% 
  mutate(term = case_match(term, 
                           "pop.speech" ~ "Populist Discourse", 
                           "pop.speech:winner1" ~ "Pop. Discourse. \nx Pres. Supporter",
                           "pop.speech:non_voter1" ~ "Pop. Discourse. \nx Non-Voter"))

# Check the work a bit -----------------------------

# Bolivia in the original lists
countries[10] 
drop_bol1 <- tidy(cong_models2_it[[10]]) %>% 
  filter(term %in% vars)

#Bolivia is indeed dropped 
cong_models2_it[[10]]$frame %>% 
  select(country) %>% 
  unique()

# Bolivia in the plotted dataframes 
drop_bol2 <- cong_models %>% 
  filter(model == "BOL", submodel == 2)

#Estimates are the same 
drop_bol1$estimate == drop_bol2$estimate

# Same check with Argentina in court models 
countries[17]
drop_arg1 <- tidy(court_models3_it[[17]]) %>% 
  filter(term %in% vars)

# Arg is dropped 
cong_models3_it[[17]]$frame %>% 
  select(country) %>% 
  unique()

# Arg in plotted dfs
drop_arg2 <- court_models %>% 
  filter(model == "ARG", submodel == 6)

# Ests are the same 
drop_arg1$estimate == drop_arg2$estimate


# Final Plots ----------------------------------------

small_multiple(cong_models) +
  geom_hline(yintercept = 0,
             colour = "red",
             linetype = 2) +
  labs(title = "Models of Support for Closing Congress", 
       x = "Country Dropped") +
  scale_color_hue(name = "Model") +
  theme_light() +
  theme(axis.text = element_text(color = "black"), 
        strip.text = element_text(color = "black"))

ggsave(here("figures/congress_iterative.pdf"), height = 6, width = 8)

#path <- paste("C:/Users/L03534594/Dropbox/Apps/Overleaf/Populism and Support for Executive Aggrandizement - CPS 2023/rr_figs/")

#ggsave(paste0(path, "si_cong_iterative.pdf"), height = 6, width = 8)


small_multiple(court_models) +
  geom_hline(yintercept = 0,
             colour = "red",
             linetype = 2) +
  labs(title = "Models of Support for Closing the Court", 
       x = "Country Dropped") +
  scale_color_hue(name = "Model") +
  theme_light() +
  theme(axis.text = element_text(color = "black"), 
        strip.text = element_text(color = "black"))

ggsave(here("figures/court_iterative.pdf"), height = 6, width = 8)

#ggsave(paste0(path, "si_court_iterative.pdf"), height = 6, width = 8)




